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Dynamical networks with time delays can pose a considerable challenge for mathematical analysis. 
Here, we extend the approach of generalized modeling to investigate the stability of large networks 
of delay-coupled delay oscillators. When the local dynamical stability of the network is plotted 
as a function of the two delays then a pattern of tongues is revealed. Exploiting a link between 
structure and dynamics, we identify conditions under which perturbations of the topology have a 
strong impact on the stability. If these critical regions are avoided the local stability of large random 
networks can be well approximated analytically. 
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Dynamical networks with time-delays (delay networks) 
have many applications in diverse range of fields from 
physics and biology [1-8] . In particular in biological sys- 
tems, both the functional forms governing individual dy- 
namical elements and the precise topology of interdepen- 
dencies are often uncertain. For making progress it is 
therefore crucial to gain a general understanding which 
properties of the local dynamics and the coupling topol- 
ogy have a strong impact on the system-level dynamics. 

The analysis of delay systems is challenging, because 
even a single delay-differential equation (DDE) consti- 
tutes an infinite-dimensional dynamical system. In the 
past, investigations of delay networks in continuous time 
have primarily focused on small or structurally simple 
systems [9-13]. Numerical explorations of larger net- 
works emphasized different behavior can be observed de- 
pending on the coupling topology [14], which can be re- 
lated to the spectrum of the graph Laplacian [15, 16]. 
Using a different approach the effect of random delays 
was analyzed successfully by a mean-field approximation 
[17]. 

For a variety of dynamical systems without delay it 
was recently shown that the dynamics of large networks 
can be analyzed efficiently by the approach of general- 
ized modeling [18, 19]. In particular this approach was 
recently applied to systems from ecology [20, 21] and cell 
biology [22, 23]. In this context, the ability to incorporate 
delays into generalized models is highly desirable. 

We consider a class of models where time delays ap- 
pear both in each network node and in the coupling of 
nodes. This is inspired by earlier works that showed that 
the interplay of two different delays in a single node [24] , 
or in the coupling [25], gives rise to rich dynamics. Both 
intra-node and coupling delays are present in ecological 
metacommunities (times needed for maturation of indi- 
viduals and migration between patches) and systems bi- 
ology (protein assembly time, active/passive transport 
times) . 

Here we use generalized modeling to obtain an expres- 
sion governing the stability of all stationary states in 
this class of system. A numerical analysis reveals a rich 



pattern of tongues of different instabilities in the space 
spanned by the two delays. We then propose an approx- 
imation, allowing for the analytical investigation of the 
pattern of instability. Thereby, we identify the condi- 
tions under which small variations in the network topol- 
ogy have a strong impact on stability. If these regions 
are avoided then the local stability of large networks of 
delay-coupled delay oscillators can be well approximated 
analytically. 



DELAY-COUPLED DELAY NETWORKS 

We consider networks of N nodes and K bidirectional 
links. The topology of these networks is captured by the 
adjacency matrix A, such that Ay is 1 if nodes i and 
j are connected and otherwise. Each node i has an 
internal dynamical variable Xi , representing for instance 
the abundance of an ecological population or the density 
of mRNA molecules. Hence, Xi changes due to internal 
dynamics in i and coupling to the neighbor's variables 
Xj according to 

X t = G{XJ) - L{Xi) + ]T Ay(F(X|) - F(Xi)), (1) 
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where r and S denote internal and travel-time delays and 
G, L and F are positive functions describing gain, loss, 
and coupling, respectively. In our analytical treatment 
we do not restrict these functions to specific functional 
forms, but consider formally the whole class of models, 
which includes several well-studied examples such as the 
Mackey-Glass [26] and the Ikeda model [27]. 

GENERAL STABILITY ANALYSIS 

The challenge that we address in the following is to 
determine the local stability of an arbitrary (positive) 
steady state Xi* . In the physical system such states 
may either correspond to stationary points (e.g. oscil- 
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lation death) or phase-locked oscillations described in a 
co-moving frame. 

Previous works on generalized models showed that it 
is possible to express the stability of arbitrary steady 
states as functions of easily interpretable parameters if 
a certain normalization is used [18, 19]. We therefore 
introduce normalized variables, Xi = Xi/X*, and nor- 
malized functions, f(Xi) — F(xiX*)/F(X*), denoted in 
lower case. Using X* = X T * = X s * , we write 

±i = a(g(xl) - l{ Xi )) + /3J2 - f(xi)), (2) 



where a = G(X*)/X* = L(X*)/X* and (3 = F(X*)/X* 
can be interpreted as parameters describing characteris- 
tic per-capita rates of growth and transport. Finally, we 
set a = 1 by normalization of the time scale. 

The stability of the normalized steady state x* = 1 
can be determined from a local linearization given by the 
Jacobian matrix J with Jy = dXi/dXj. Close to the 
steady state, small perturbations can be decomposed into 
eigenvectors vi of the Jacobian [28]. The time evolution 
of a perturbation y along the eigenvector vi , then follows 
a locally exponential trajectory such that 



dy 



= e" AlT , 



(3) 



where A/ is an eigenvalue of the Jacobian, corresponding 
to the eigenvector «/. Using this relation we capture the 
response to a given perturbation by a Jacobian with 
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where dashes denote derivatives with respect to the argu- 
ment in x* and A is a self-consistent eigenvalue of J (A). 

The steady state is stable if all of these eigenvalues 
have negative real parts. Conversely, the steady state is 
unstable if at least one such eigenvalue has positive real 
part, such that there is a perturbation that grows in time. 

We emphasize that Eq. (4) states the Jacobian of 
all steady states in all models of the form of Eq. (1). 
The Jacobian is written as a function of scalar quan- 
tities that can be interpreted as unknown parameters. 
Specifically, the parameters <?', V , and /' are logarithmic 
derivatives of the original functions (e.g. /' = df /dx\ x = 
d\og(F) / d\og(X)\ Sf ) , which are known as elasticities and 
are used in many fields because they often have an in- 
tuitive interpretation in the context of the application 
[29]. 

Because the present analysis focuses primarily on the 
effect of topology, we restrict ourselves to the case g' = 
— 1, V = 0, /' = 1, j9 = 1. A detailed analysis of the 
effect of parameters will be published separately. 

In DDEs the computation of eigenvalues is complicated 
by the explicit appearance of the eigenvalue A in J, which 
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FIG. 1. (Color online) Computation of the number of eigen- 
values with positive real part. Left: Location of the largest 
eigenvalues of an example Jacobian in the complex plane (cir- 
cles). The algorithm follows a rectangular contour C (yel- 
low), taking discrete steps (dots). The size of C, marked o, b, 
is chosen sufficiently large that all positive eigenvalues must 
lie within contour (Gershgorin's theorem). Right: The num- 
ber of positive eigenvalues is found as the winding number of 
arg P(X) as C is followed in positive (counter clockwise) direc- 
tion. Colors (black/red/blue/gray) make different segments 
of C. 



turns the characteristic polynomial P(X) into an implicit 
transcendental equation. We therefore follow the ap- 
proach of [30] and test for eigenvalues with positive real 
parts (EVPs) using Cauchy's argument principle: For an- 
alytic functions P(X), the number of roots inside a con- 
tour C is Nc = T^AcargP(A) where AcargP(A) is the 
winding number of P(X) on C (Fig. 1). The total num- 
ber of EVPs can therefore be computed by applying the 
Cauchy principle to a contour in the positive half-plane 
that is chosen so large that all EVPs must lie within the 
contour. In the present Letter we use a rectangular con- 
tour covering the interval [0, a] in the real and [— b, b] in 
the imaginary direction, where a and b are estimated us- 
ing Gershgorin's theorem [31]. For efficient computation 
of the winding number we use an adaptive step-size algo- 
rithm counting the crossings of odd multiples of it. We 
iterate along C by stepping from a point Xi to the sub- 
sequent point Ai+i = Xi + hv, where v is a unit vector 
along C and h is the current step size. The step is ac- 
cepted if a) A := D(axgP(zi) ~ ar gf (^i+i))) < e = 0.1 
and b) L>(argP(z 4 ) - arg P((z; + z i+ i)/2)) < A where 
D(x) = \x\ mod 2-7T. Otherwise, the step is rejected and 
h — s- h/2. After each successful step h — > /imax(2, e/A). 



STABILITY IN FULLY-CONNECTED 
NETWORKS 

We now explore the effect of topology and delay on lo- 
cal stability of steady states. For understanding the effect 
of delays on dynamical stability we first consider a fully- 
connected network. The stability of this network is ex- 
plored numerically by generating an ensemble of param- 
eter sets, where the delays r and 5 are drawn randomly 
from a uniform distribution. The stability of the steady 
states corresponding to the sample parameter sets is then 
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FIG. 2. (Color online) Local stability and bifurcations of networks. The parameter space in the panels is sampled uniformly. 
Stable states are marked in green, whereas unstable states are not shown. Panels show stability in the fully-connected network 
with N = 10 nodes and K = 45 links (a), random trees with N = 10 and K = 9(b), and random networks N = 10 and K = 15 
(c), and ./V = 100 and K = 500 (d). Areas of instability (white) are sharply delineated if a single topology is considered (a). In 
network ensembles fuzzy regions appear, where stability is topology dependent (b,c). In case of random trees these regions are 
bounded by the bifurcation lines of star (dashed) and chain (solid) topologies (c). However, even in large random networks with 
many links the regions of instability have relatively sharp boundaries, which can be approximated by analytical bifurcation 
lines for the fully-connected network of identical degree d = 10 (d, solid). 



evaluated by the method described above. Because of the 
efficiency of this method the evaluation of the ensemble 
only requires minutes of computational time. 

The result of the sampling analysis reveals a pattern 
of tongues of instability, which are separated by stable 
channels located around the resonant delays r = nS with 
integer n (Fig. 2a). 

We note that these results differ from those found in 
systems with two delays of the same type [24, 25], where 
resonance takes the form mr = nS. 

Qualitative changes in the dynamics, including changes 
in the stability can only occur when the system undergoes 
a bifurcation. For the fully-connected network we com- 
puted the local bifurcation points explicitly by numeri- 
cal continuation of the bifurcation condition Re(A) = 0. 
Some of the bifurcations that are thus revealed mark 
changes of stability on the edge of the tongue, whereas 
others correspond to qualitative transitions within the 
unstable region. 

Considering the tongues in Fig. 2a further, one no- 
tices that the tips of the tongues are located on the ver- 
tices of a square lattice. A similar symmetry was already 
observed previously in simpler models [32] and can be 
explained as follows: Because the edge of a tongue is 
a local bifurcation, the Jacobian has to have a purely 
imaginary eigenvalue A = iu>. Moreover, a delay param- 
eter, say r, can only appear in the Jacobian in factors 
of the form e~ Ar . Therefore, increasing the delay by a 
multiple of 2tt/uj leaves the Jacobian invariant because 
e -iw(r+27r/a;) _ g-wr^ j e ^ when starting from a bifurca- 
tion point (r, 5, ui), increasing either of the two delays by 
27r/w must lead to another bifurcation point. This proves 
that not only the tips of the tongues, but every bifurca- 
tion point is part of a square lattice. For understanding 
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FIG. 3. (Color online) Lattice symmetry in the bifurcation 
points. Top: Blowup of a part of the bifurcation diagram, 
Fig. 2a. Bottom: Selected eigenvalues of the Jacobian (differ- 
ent colors) along a one dimensional cut at 8 — 2 (black line, 
top). Shaded regions mark values of r where the system with 
S — 2 is unstable. 



why the periodicity is most visible in the tips, consider 
that different points on the edge of a single tongue gen- 
erally differ in uj and therefore also in the correspond- 
ing lattice constant. Therefore, the tongues as such do 
not reappear identically at higher delays but become dis- 
torted. We emphasize that this symmetry cannot be ex- 
tended to points in parameter space that are not bifur- 
cation points (see Fig. 3). 
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STABILITY IN RANDOM TREES 

We investigate the impact of network topology on the 
stability first by considering an ensemble of random tree 
topologies. For any given tree the analysis described 
above reveals sharply delineated tongues of stability. The 
number and location of these tongues depends on the spe- 
cific topology under consideration. For visualizing the 
differences and similarities between tongue patterns, we 
repeat the sampling analysis described above but also 
draw a random tree topology for every individual sample 
point. Because results from many different topologies are 
thus superimposed in Fig. 2b, the edges of the tongues 
become fuzzy. The fuzzy areas thus mark parameter re- 
gions where the stability depends on the network topol- 
ogy. Nevertheless, large region exist in which the stability 
remains independent of topology. 

We found empirically that the stability boundaries for 
all tree topologies fall between those for the linear chain 
and the star topology, which are the least and most dy- 
namically stable topologies, respectively. In some places 
the stability boundaries of even these extreme cases co- 
incide so that the corresponding transitions occur inde- 
pendently of the specific tree topology. 

In contrast to the fully-connected network, the peri- 
odicity in the travel-time delay, 5, is reduced to tt/uj for 
trees. The resulting additional tongues cover the stable 
channel around r = S, present in the fully-connected net- 
work. 



STABILITY IN RANDOM NETWORKS 

For further exploration of the effect of topology and 
understanding the appearance of additional tongues of 
instability in the random trees, we now consider an en- 
semble of Erdos-Renyf random graphs, with fixed number 
of links, L. 

Below, we distinguish two 27r-periodic sets of tongues: 
the off-diagonal set (OS), which is already present in 
the fully-connected network, and the diagonal set (DS), 
which appeared in the trees. By visual inspection of the 
phase diagrams for different topologies we observed that 
the DS is present in some networks, but is absent in oth- 
ers. 

While the 2n periodicity within the sets is guaranteed 
by analytical arguments, the relative offset between the 
sets can depend on topology. The 7r-periodicity, observed 
in the random trees, requires that the DS is shifted rela- 
tive to the OS exactly by (0, tt/ui) in (r, 5). We find this 
particular offset whenever the network topology is bipar- 
tite, i.e. when the network can be colored with two colors 
such that no link connects nodes of the same color. The 
observed 7r-periodicity therefore arises because all trees 
meet the condition of bipartiteness. 



The dependence of the stability on the topology noted 
earlier has to stem from differences either in the num- 
ber or location of the tongues of stability. Given two 
network topologies, parameter regions with different sta- 
bility properties appear because a) the tongues found 
for one of the topologies are shifted with respect to the 
tongues found for the other topology, or b) tongues of 
the DS are present in one of the topologies but absent 
in the other. In comparison, the effect of a) is relatively 
minor, causing for instance the small regions of topol- 
ogy dependence in Fig. 2b. By contrast, if K is tuned 
to a value where the DS appears, stability can be topol- 
ogy dependent in large regions of the parameter space 
(Fig. 2c). 

Even when ensembles of large random networks are 
considered, we observe that the tongues of instability re- 
main relatively sharply delineated (Fig. 2d) unless pa- 
rameters are tuned to regions of parameter space where 
new tongues appear (Fig. 2c). This implies that in the 
ensemble considered here, the dynamical stability is to a 
large extend independent of the specific network topol- 
ogy. We note that the stability of large and dense random 
networks closely matches the results for a fully-connected 
network with the same mean degree (Fig. 2d). For the 
special case of degree-homogeneous networks, this obser- 
vation is explained below. 

ANALYTICAL THEORY 

For gaining an analytical understanding of the re- 
sults presented above, we consider the case of degree- 
homogeneous networks, in which all nodes have the same 
number of connections. For these networks the eigenval- 
ues of the Jacobian from Eq. (4) are given by the implicit 
equation 

A = J d (A) + cJ°(A), (5) 

where c can be any eigenvalue of the adjacency matrix, 
which we denote as topological eigenvalues. Bifurcations 
are characterized by the presence of a purely imaginary 
eigenvalue iu). This leads to 

- g' cos(0) -I'- d/3f + cfif cob(V), (g) 
w = — g' sm{4>) — cj3f sin(?/>). 

where we used Eq. (4) and <j> := ljt, tJj :— ojS. 

For analyzing the bifurcation condition we enumer- 
ate its solution branches. First, we note that for 
each solution triplet ((f>,ip,u), there exist another solu- 
tion (—0, —ip, —ui) corresponding to the identical tongue. 
Therefore, we only consider solutions with ui > 0. Sec- 
ond, given a solution {(p, ip, lu) other solutions are found 
at (0 + 2nr, ip + 2tts, ui), where r, s are integers enumer- 
ating the tongues. Finally, if there is a triplet (</>, ip,u>) 
that solves Eq. (6) for c < 0, there must be a triplet 
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(<f>, if) + 7T, ui) that solves modified equations in which c is 
replaced by |c|. We can therefore enumerate the solutions 
with negative c by half-integer values of s in the modified 
system. 

Considering the first tongue (r = s = 0) we find 



i> = 27r±cos~ 1 (p(</>)) 

w = -g' sin(0) T M/3/Vl-pW, 



where 



p(4>) 



V — g' cos(0) 



(7) 



(8) 



and the codomain of cos -1 is [0, 7r]. These equations 
provide a representation of the tongues depending on <f>. 

We define the tip of a tongue as the bifurcation point 
recurring with the smallest lattice constant. The tips are 
thus characterized by maximal values of w, such that 



Ttip 



cos 1 (q) + 2rn 



(\c\Pf- 
cos _1 (g) 



(2s + 1)tt 



(9) 



with g = + Z')/(|c|/3/' - ff')- Eq. (9) shows that 

half-integer values of s, or equivalently topological eigen- 
values c < 0, correspond to the DS, while integer values of 
s correspond to the OS. Further, Eq. (9) requires |g| < 1 
and hence 



Id > d + 



g' + 1 



(10) 



Topological eigenvalues c violating this condition cannot 
satisfy the bifurcation condition and hence do not corre- 
spond to tongues of instability. The positive c creating 
the OS, and the negative c creating the DS are there- 
fore separated by a forbidden region in which topological 
eigenvalues do not create tongues. If a change of param- 
eters causes a c to enter this region, the corresponding 
set of tongues vanishes as the corresponding bifurcation 
lines shift to infinite delays. 

For understanding the appearance of 7r-periodicity in 5 
observed in bipartite networks, consider that all bipartite 
networks have symmetric topological spectra [33] . There- 
fore, for every c > there is a symmetric c < 0, such that 
the 7r-periodicity appears according to Eq. (6). 

Above we observed that the stability of large and dense 
random networks closely matches the results for the fully- 
connected network with the same mean degree. This de- 
pendence of the similarity on the mean degree can be 
explained in degree-homogeneous networks as the largest 
eigenvalue in these networks is c max = d. The OSs cre- 
ated by the leading eigenvalues of two different degree 
homogeneous networks therefore match if the networks 
have the same degree. Further, the observed similarity 
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FIG. 4. (Color online) Comparison between simulations and 
stability analysis for a ring of 10 Mackey-Glass systems, a) 
Bifurcation lines of the topological eigenvalues 2 (solid black) 
and —2 (dotted red), (b-d) Simulation results for two neigh- 
boring nodes with r = 0.28,5 = 0.6 (b), r = 0.6,5 = 0.6 (c) 
and r = 1.8,5 = 0.9 (d). Delay values inside tongues of the 
eigenvalue 2 give rise to in-phase oscillations, values inside 
tongues of eigenvalue —2 give rise to anti-phase oscillations. 
Simulations are performed with the pydelay-package [34]. 



requires that no other sets of tongues are present, which 
seems to be generally true for sufficiently large and dense 
random graphs and fully-connected networks. This ex- 
plains that fully-connected networks with appropriately 
chosen mean degree offer a good approximation for the 
stability in a large class of random graphs with suffi- 
ciently narrow degree distribution. 



NUMERICAL SIMULATIONS 

For confirming the results of the generalized model and 
illustration of the dynamical implications of bifurcation 
lines we simulate 10 Mackey-Glass systems coupled in a 
ring topology. The simulations use the equations 



Xi = 



aXJ 



i + (x;y 



-cX i +e(Xf_ 1 + Xf +1 -2X i ), (11) 



with a = 2, b = 10, c = 1, e = 10, corresponding to the 
general parameters g' = —4,1' = l,f = l,/3 = 10. The 
two topological eigenvalues c\ = 2 and c\q — —2 satis- 
fying Eq. 10 give rise to one DS and one OS of tongues. 
The dynamical implications of the tongues is illustrated 
in Fig. 4. Choosing parameters inside the OS results in 
in-phase synchronized oscillations, parameters inside the 
DS, result in anti-phase oscillations, while parameter val- 
ues outside the tongues result in stationary dynamics. 
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DISCUSSION 

In this paper we studied the local dynamics of a large 
class of networks of delay-coupled delay differential equa- 
tions. For these we showed that the stability of steady 
states is governed by two distinct sets of tongues of insta- 
bility. Further, exploiting a link between network struc- 
ture and dynamics, we derived an analytical expression 
for the tongues in degree-homogeneous nets. 

Our results were obtained using the approach of 
generalized modeling that constitutes a local analysis. 
Presently extensions of this approach to nonlocal dynam- 
ics are being developed [35]. However, even the present 
local analysis can reveal some insights in nonlocal phe- 
nomena. For instance it is known that chaotic dynamics 
are generically present close to double Hopf bifurcations 
that are found at the intersection of tongues [36, 37]. 

The main result of our exploration is that the stability 
in large random networks is relatively independent of the 
network topology. This result seems to be in conflict 
with previous analysis of other systems which suggest 
that topological properties such as clustering and cliques 
can have a strong impact on the dynamics [14, 16]. The 
apparent disagreement is resolved if one considers that 
such topological properties are exceedingly rare in the 
ensemble of random networks considered here. 

Beyond the numerical results for random networks, our 
analytical calculations showed that the dynamics are in 
general only strongly dependent on topology if a variation 
of the topology brings new tongues of instability into ex- 
istence. This is the case if the variation shifts eigenvalues 
of the adjacency matrix into a certain range. This obser- 
vation is consistent with previous results as topological 
properties that have been implicated as influencing the 
dynamics are closely linked to spectral properties of the 
networks [38-41]. We can therefore conjecture that the 
strong topology dependence of the dynamics observed 
previously is linked to specific eigenvalues of the adja- 
cency matrix that appear for instance in networks con- 
taining many clusters or cliques. It was recently shown 
that the appearance of these eigenvalues is not linked to 
the clusters or cliques themselves, but rather to specific 
local symmetries that typically accompany them [42] . 

If true, the conjecture stated above opens an intriguing 
possibility: At least in the class of networks considered 
here, the local dynamics of very different networks with 
the identical mean degree may indeed be very similar 
except for additional instabilities which are caused by 
local symmetries within the networks. 
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